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Abstract 

This  paper  addresses  the  design  and  computation  of  a  guidance  law  for  a  transfer  mission  from  an  orbit  near  the  Earth  to  a  halo  orbit 
around  the  libration  point  L2  in  the  Sun-Earth  system.  The  guidance  law,  which  is  designed  based  on  receding  horizon  control  and  com¬ 
pensates  for  launch  velocity  errors  that  are  introduced  by  inaccuracies  of  the  launch  vehicle,  is  solved  using  the  generating  function 
method.  During  the  design  of  the  closed-loop  guidance  law,  the  entire  transfer  mission,  which  is  considered  a  nonlinear  optimal  control 
problem,  is  evaluated  to  obtain  a  nominal  reference  trajectory.  Using  the  launch  velocity  errors  and  the  uncertainty  of  the  model,  a  space¬ 
craft  controlled  by  the  proposed  guidance  law  tracks  the  reference  trajectory.  Furthermore,  the  original  Riccati  differential  equation  in 
the  receding  horizon  control  algorithm  is  replaced  by  an  equivalent  convenient  form  of  the  Riccati  differential  equation  that  is  based  on 
the  generating  function.  The  high-efficiency  solution  of  the  equivalent  equation  avoids  the  online  direct  integration  of  the  original  Riccati 
differential  equation,  which  significantly  increases  the  computational  efficiency  for  the  receding  horizon  control  problem.  Numerical  sim¬ 
ulations  using  a  nonlinear  bicircular  four-body  model  demonstrate  the  capabilities  of  the  proposed  receding  horizon  guidance  law  for  the 
transfer  mission.  In  addition,  the  generating  function  method  improves  the  computational  efficiency  by  at  least  one  order  of  magnitude 
over  the  backward  sweep  method  in  solving  the  receding  horizon  control  problem. 
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1.  Introduction 

The  exploration  of  the  libration  points  in  space  has  been 
accompanied  by  a  rich  and  growing  literature  on  the 
design,  modeling,  and  control  of  a  variety  of  orbits  (Gomez 
et  al.,  2003;  Farquhar  et  al.,  2004;  Wang  et  al.,  2007).  The 
halo  orbit,  which  is  a  special  three-dimensional  periodic 
libration  point  orbit,  is  useful  for  performing  experiments 
and  observations  about  the  solar  and  space  structure  and 
providing  continuous  communications  between  the  Earth 
and  the  far  side  of  the  Moon  (Farquhar,  1970a, b;  Break- 

*  Corresponding  author.  Tel./fax:  +86  411  84706574. 

E-mail  addresses:  hjpeng@dlut.edu.cn  (H.  Peng),  qgao@dlut.edu.cn 
(Q.  Gao),  wuzhg@dlut.edu.cn  (Z.  Wu),  zwoffice@dlut.edu.cn  (W.  Zhong). 


well  et  al.,  1974;  Farquhar  et  al.,  1980;  Howell  and  Perni- 
cka,  1993;  Romagnoli  and  Circi,  2010;  Zanzottera  et  al., 
2012).  Several  space  missions  have  involved  libration 
points,  including  ISEE-3,  Wind,  SOHO,  ACE  and  Genesis 
(Xu  and  Xu,  2008).  The  scientific  goals  of  these  missions 
were  to  examine  the  structure  of  the  solar  wind  near  the 
Earth  and  the  shock  wave  that  forms  the  interface  between 
the  solar  wind  and  the  Earth’s  magnetosphere,  to  monitor 
the  solar  wind  and  associated  phenomena  near  the  libra¬ 
tion  point,  and  to  study  the  energetic  particles  from  the 
solar  wind,  the  interplanetary  medium,  and  other  sources 
(Farquhar,  1970a, b;  Breakwell  et  al.,  1974;  Howell  and 
Pernicka,  1993;  Romagnoli  and  Circi,  2010;  Zanzottera  et 
al.,  2012).  Therefore,  the  problem  of  transferring  a  space¬ 
craft  from  the  Earth  to  libration  point  orbits,  especially 
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halo  orbits,  is  an  important  problem  in  the  field  of  explora¬ 
tion  of  the  libration  point  space  environment. 

Because  of  the  fundamental  advantages  of  the  space 
environment  in  libration  orbits,  transfers  to  these  libration 
point  orbits  have  been  widely  studied  (Gomez  et  al.,  2005; 
Tantardini  et  al.,  2010).  Most  of  the  previous  studies  on 
transfer  missions  have  focused  on  dynamical  system  mod¬ 
els  and  the  type  of  thrust  needed  for  the  transfer.  The  cir¬ 
cular  restricted  three-body  problem  (CRTBP)  model  is  the 
dynamical  model  that  is  most  commonly  used  as  a  founda¬ 
tion  for  designing  transfer  orbits  (Gomez  et  al.,  2005;  Tan¬ 
tardini  et  al.,  2010;  Li  and  Zheng,  2010;  Serban  et  al., 
2002).  Conley  (1968)  studied  transit,  non-transit  and 
asymptotic  orbits  around  the  Lagrange  point  based  on 
the  CRTBP  model.  The  effects  of  the  Sun  in  the  Earth- 
Moon  system  are  ignored,  and  the  effects  of  the  Moon  in 
the  Sun-Earth  system  are  also  ignored.  Hiday-Johnston 
and  Howell  (1994)  formulated  a  strategy  to  design  optimal 
time-fixed  impulsive  transfers  between  three-dimensional 
libration-point  orbits  based  on  the  elliptic  restricted 
three-body  model.  The  four-body  dynamic  model  was 
developed  to  consider  the  effects  of  the  Sun  or  the  Moon 
in  the  design  of  transfer  orbits  (Salmani  and  Biiskens, 
2011;  Cabette  and  Prado,  2008;  Assadian  and  Pourtakd- 
oust,  2010).  The  Earth-Moon-Sun-Satellite  system  is  used 
to  study  the  transit  trajectory  properties  in  restricted  three- 
and  four-body  problems  (Circi,  2012). 

This  paper  uses  the  so-called  bicircular  four-body 
model,  which  places  the  Sun  and  Earth  in  circular  orbits 
about  their  barycenter  and  places  the  Moon  in  a  circular 
orbit  about  the  barycenter  of  the  entire  system  (the  four- 
body  dynamical  model  was  developed  by  Salmani  and  Bus- 
kens,  2011).  The  advantage  of  this  model  is  that  it  is  more 
realistic  than  both  the  CRTBP  of  the  Sun-Earth  system 
and  the  CRTBP  of  the  Earth-Moon  system.  The  main  dis¬ 
advantage  of  this  model  is  that  it  destroys  the  periodic 
solutions  that  exist  in  the  CRTBP  model  (Salmani  and  Bus- 
kens,  2011).  In  addition,  the  bicircular  model  is  time-peri¬ 
odic  with  two  natural  frequencies,  and  the  periodic  orbits 
are  replaced  by  quasi-periodic  orbits.  In  this  paper,  the 
halo  orbit  around  the  libration  point  L2  in  the  Sun-Earth 
system  is  considered  the  objective  orbit,  so  the  CRTBP 
model  is  considered  to  generate  the  periodic  halo  orbit  to 
overcome  the  lack  of  periodic  solutions  described  above. 

The  choice  of  a  dynamical  model  is  dependent  on  the 
requirements  of  the  research  mission.  A  more  complex 
N-body  numerical  solution  and  real  JPL  ephemeris  data 
have  higher  accuracies  than  the  CRTBP  model  and  the 
bicircular  four-body  model  for  orbit  design.  Because 
the  purpose  of  this  paper  is  to  design  and  compute  an 
online  guidance  law  for  a  transfer  mission,  the  more  com¬ 
plex  dynamical  model  requires  additional  online  computa¬ 
tional  time.  In  addition,  solar  pressure,  initial  velocity 
errors  and  the  inaccuracies  of  the  thrusts  will  greatly  influ¬ 
ence  the  dynamical  model.  Therefore,  the  appropriate 
accuracy  of  the  dynamical  model  is  used  in  this  paper, 
and  the  influences  of  the  modeling  uncertainty  will  be 


repressed  by  the  online  guidance  law  and  not  by  the  more 
complex  dynamical  model. 

Another  important  problem  is  the  type  of  thrust  that  is 
used  for  the  transfer  to  the  libration  point  orbit.  Indeed, 
the  dynamical  systems  approach  (Mingo tti  et  al.,  2012, 
2007;  Kechichian,  2001)  has  been  employed  as  a  low- 
energy  transfer  method  to  transfer  a  spacecraft  to  a  halo 
orbit  with  its  stable  manifold.  The  Genesis  trajectory  was 
the  first  mission  to  be  fully  designed  using  dynamical  sys¬ 
tem  theory  (Serban  et  al.,  2002;  Howell  et  al.,  1997).  How¬ 
ever,  only  the  stable  manifold  associated  with  the  large 
halo  orbit  is  connected  with  the  parking  orbit  near  the 
Earth.  Furthermore,  a  trajectory  correction  maneuver 
should  be  implemented  during  the  long-duration  transfer 
to  halo  orbit.  The  two  most  important  approaches  for  tra¬ 
jectory  correction  are  the  impulsive  maneuver  (Li  and 
Zheng,  2010;  Rausch,  2005)  and  low- thrust  propulsion 
(Ozimek  and  Howell,  2010;  Dellnitz  et  al.,  2009).  In  gen¬ 
eral,  the  impulsive  maneuver  uses  a  parameter  optimization 
problem  to  design  the  thrust  and  direction  variables.  Low- 
thrust  propulsion  for  transfer  to  a  halo  orbit  can  be  mod¬ 
eled  as  a  continuous  optimal  control  problem  and  is  more 
difficult  than  the  parameter  optimization  problem.  While 
high-thrust  impulsive  systems  have  not  received  sufficient 
attention  over  the  last  decade,  many  studies  have  addressed 
low-thrust  spacecraft  trajectory  optimization.  This  paper 
models  the  transfer  problem  based  on  the  continuous 
low-thrust  as  an  optimal  control  problem  that  is  solved 
using  a  symplectic  approach  (Peng  et  al.,  2011). 

In  practice,  planning  the  optimal  transfer  trajectory  is 
not  sufficient  to  implement  a  transfer  mission.  Deviations 
from  the  required  thrust  are  introduced  by  inaccuracies 
in  the  engine  and  must  be  corrected  by  subsequent  maneu¬ 
vers.  In  addition,  the  transfer  mission  to  the  libration  point 
orbit  is  extremely  sensitive  to  launch  errors.  Therefore,  the 
guidance  strategy  should  compensate  for  all  of  the  distur¬ 
bances  that  cause  instabilities  in  the  spacecraft’s  trajectory. 
Guidance  methods  such  as  predictor-corrector  guidance 
and  reference  trajectory  guidance  can  be  used  to  correct 
deviations  from  the  reference  nominal  trajectory  (Tian 
and  Zong,  2011;  Arrieta-Camacho  and  Biegler,  2005). 
The  predictor-corrector  guidance  method  can  handle  large 
dispersions  and  accommodate  severe  off-nominal  condi¬ 
tions.  However,  for  highly  constrained  guidance  scenarios, 
the  computational  requirements  and  convergence  guaran¬ 
tees  are  the  most  difficult  problems  (Tian  and  Zong, 
2011).  The  reference  trajectory  guidance  method  is  used 
to  track  a  reference  trajectory  and  is  implemented  with 
an  online  closed-loop  guidance  law.  Receding  horizon  con¬ 
trol  has  been  used  as  a  reference  trajectory  guidance 
method  for  precision  entry  guidance  missions  (Lu,  2000, 
1999).  When  using  receding  horizon  control,  the  optimal 
control  problem  is  solved  over  a  shorter  moving  horizon, 
and  the  Riccati  differential  equation  needs  to  be  integrated 
backward  over  this  finite  horizon  for  every  time  point 
(Kwon  and  Pearson,  1977).  However,  the  online  computa¬ 
tional  burden  associated  with  solving  the  Riccati  differen- 
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tial  equation  is  an  impediment  for  practical  applications  of 
the  method  (Lu,  2000).  Nonetheless,  the  methods  described 
above  are  important  and  serve  as  the  foundation  for  the 
development  of  practical  and  feasible  control  methods. 

Because  of  the  issues  described  above,  a  highly  efficient 
numerical  algorithm  to  reduce  the  online  integration  bur¬ 
den  of  the  Riccati  differential  equation  has  become  a  key 
problem  for  the  real-time  implementation  of  receding  hori¬ 
zon  control.  A  typical  method  that  is  based  on  the  approx¬ 
imate  state  and  control  at  discrete  time  points  (Lu,  2000) 
avoids  the  online  integration  of  the  Riccati  differential 
equation.  However,  the  numerical  accuracy  of  the  receding 
horizon  control  problem  is  dependent  on  the  approximate 
state  and  control  variables.  Because  there  are  no  approxi¬ 
mations  for  the  state  and  control  variables,  the  advantage 
of  the  online  integration  of  the  Riccati  differential  equation 
is  its  high  accuracy.  Therefore,  this  paper  first  avoids  the 
direct  integration  of  the  original  Riccati  differential  equa¬ 
tion  and  then  provides  an  equivalent  equation  that  is 
expressed  by  the  generating  function  (Wu  and  Zhong, 
2009;  Wu  and  Mesbahi,  2012).  A  highly  efficient  numerical 
method  is  proposed  for  solving  the  equivalent  equation. 
Not  only  can  the  numerical  accuracy  for  receding  horizon 
control  be  guaranteed,  but  the  online  computational  effi¬ 
ciency  can  also  be  greatly  increased.  The  underlying  reason 
is  that  in  the  detailed  implementation  of  the  proposed 
numerical  algorithm,  the  solution  to  the  equivalent  equa¬ 
tion  can  be  obtained  by  the  combination  and  separation 
of  formulas  of  the  generating  function.  This  numerical 
algorithm  avoids  the  online  direct  integration  of  the  origi¬ 
nal  Riccati  differential  equation.  Meanwhile,  the  proposed 
numerical  algorithm  can  ensure  the  closed-loop  stability  of 
the  receding  horizon  control  method. 

This  paper  is  organized  as  follows.  After  a  brief  intro¬ 
duction  to  the  bicircular  four-body  problem  of  the  Sun- 
Earth  system,  the  objective  periodic  halo  orbit  is  computed 
using  the  differential  correction  algorithm  in  Section  2.  In 
Section  3,  the  optimal  control  problem  of  the  low-thrust 
transfer  mission  to  the  halo  orbit  is  established  and  solved 
using  a  symplectic  algorithm.  In  Section  4,  closed-loop 
guidance  based  on  the  receding  horizon  control  method 
is  proposed  to  correct  the  trajectory  using  the  linearization 
equation  for  the  optimal  reference  trajectory  and  is  solved 
using  the  generating  function  method.  In  Section  5,  numer¬ 
ical  simulations  are  performed  with  the  nonlinear  bicircular 
four-body  model.  The  effects  of  the  initial  errors,  modeling 
uncertainty  and  the  guidance  law  parameters  are  discussed. 
The  numerical  results  demonstrate  the  validity  and  applica¬ 
bility  of  the  method.  Finally,  the  conclusions  are  given  in 
Section  6. 

2.  Dynamical  model 

The  main  problem  addressed  in  this  paper  is  the  design 
of  optimal  guidance  laws  to  transfer  a  spacecraft  from  an 
orbit  near  the  Earth  to  a  periodic  halo  orbit  around  the 
libration  point  L2  in  the  Sun-Earth  system.  To  increase 


the  precision  of  the  reference  transfer  trajectories,  the  bicir¬ 
cular  four-body  model  is  used  to  generate  the  optimal  low- 
thrust  transfer  trajectories  using  the  optimal  control 
method.  Because  the  target  orbit  is  a  libration  point  orbit, 
such  as  a  periodic  halo  orbit,  a  brief  explanation  of  the 
halo  orbit  is  given  in  this  section. 


2.1.  Bicircular  four-body  model 


The  bicircular  four-body  model  is  more  realistic  than 
both  the  Sun-Earth/Moon  and  the  Earth-Moon  three- 
body  models  (Salmani  and  Buskens,  2011).  It  describes 
the  dynamics  observed  in  systems  such  as  the  environment 
near  the  Earth  and  the  Moon  and  describes  the  motion  of  a 
spacecraft  of  negligible  mass  under  the  gravitational  force 
of  the  Sun,  Earth  and  Moon.  The  Sun  and  Earth  revolve 
in  circular  orbits  around  their  center  of  mass  (barycenter), 
and  the  Moon  moves  in  a  circular  orbit  around  the  center 
of  the  Earth  (Gomez  et  al.,  2003).  Based  on  the  assump¬ 
tions  of  the  bicircular  four-body  model,  a  rotating  refer¬ 
ence  frame  ( o,x,y,z )  is  defined  with  its  origin  at  the  L2 
point  of  the  Sun-Earth  system  such  that  the  x  unit  vector 
is  directed  from  the  Sun  towards  the  Earth,  the  y  unit  vec¬ 
tor  is  defined  as  being  normal  to  the  x  vector  in  the  plane  of 
the  primary  orbit  and  along  the  prograde  rotation  direction 
and  the  z  unit  vector  completes  the  right-handed  frame  and 
is  thus  normal  to  the  plane  of  the  primaries’  orbits.  The 
motion  equations  in  the  bicircular  four-body  model  can 
be  written  in  dimensionless  form  as  follows  (Gomez  et 
al.,  2003): 


x  =  2y  +  x-{\-n) 


jx+n)  (x  -  1  +  n)  {x  ~  xM) 


rL 


j>=  -2x  +  y- (1  -n)*-ii*-mM 

rPS  rPE 

z=  -(1  -At)^--^^--  rnM^- 


|3 


(1) 

(2) 

(3) 


where 


4s  =  (x  +  nf  +  /  +  z1  (4) 

rpE  =  (x  - 1  +  b)2  +  y2  +  z2  (5) 

fpM  =  {x-xMf+  (y-yuf+z2  (6) 


The  dot  represents  the  time  derivative  in  the  rotating 
frame,  /r  is  the  ratio  of  the  Earth’s  mass  to  the  sum  of 
the  masses  of  the  Earth  and  the  Sun,  rPS,  rPE  and  rPM  are 
the  distances  from  the  spacecraft  to  the  Sun,  the  Earth 
and  the  Moon,  respectively,  and  mM  is  the  parameter  of 
the  Moon  in  the  same  nondimensional  mass  units. 

Furthermore,  the  influence  of  the  Moon’s  motion  can  be 
expressed  by 

xM  =  Rm  cos  (coMt  +  60)  (7) 

yM  =  Rm  sin (coMt  +  0O)  (8) 
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where  the  Sun  and  the  Earth  hold  and  the  Moon  moves 
about  the  barycenter  of  the  Sun-Earth  system  on  a  circle 
with  radius  Rm  in  nondimensional  distance  units  (Gomez 
et  al.,  2003),  coM  is  the  angular  velocity  of  the  Moon  and 
60  is  the  initial  phase  of  the  Moon. 

2.2.  Halo  orbits 


tudes  of  the  orbits  range  from  1  x  105  km  to  5  x  105  km  in 
the  z  direction.  Although  the  sizes  of  these  halo  orbits  differ 
widely,  the  orbital  periods  are  nearly  the  same  (approxi¬ 
mately  180  days).  Meanwhile,  the  relative  velocities  of  the 
two  halo  orbits  are  smaller  than  in  the  case  of  two  body 
dynamics. 


The  advantage  of  the  bicircular  four-body  model  is  that 
it  gives  a  better  approximation  of  the  environment  near  the 
Earth;  thus,  a  transfer  trajectory  that  is  based  on  this 
model  will  consume  less  fuel.  However,  the  main  disadvan¬ 
tage  of  this  model  is  that  it  destroys  the  particular  and  peri¬ 
odic  solutions  of  the  CRTBP  model  (Gomez  et  al.,  2003). 
Therefore,  in  this  paper  the  CRTBP  model  is  used  to  com¬ 
pute  the  periodic  libration  orbits,  i.e.,  the  target  halo  orbits. 
The  CRTBP  model  can  be  obtained  easily  from  the  bicircu¬ 
lar  four-body  model  as  follows  (Schaub  and  Junkins,  2002) 


x  =  2 y  +  i 
y  =  -  2x3 


-(i 

rPS  rPE 

rPS  rPE 


(9) 

(10) 


3.  Optimal  low-thrust  transfer  to  libration  point  orbits 

In  this  section,  we  apply  the  reference  trajectory  guid¬ 
ance  method  to  establish  the  controlled  nonlinear  transfer 
formulation  and  then  generate  the  reference  optimal  low- 
thrust  transfer  trajectory  using  the  nonlinear  optimal  con¬ 
trol  method. 


3.1.  Problem  formulation 

The  transfer  mission  from  the  near-Earth  orbit  to  the 
halo  orbit  is  formulated  as  a  nonlinear  optimal  control 
problem.  The  objective  of  the  transfer  mission  is  to  mini¬ 
mize  the  energy  consumption  for  a  fixed  flight  time  (Ross, 
2006),  i.e., 


z=-(  i-ri-L-n*  (li) 

'PS  'PE 

Halo  orbits  are  special  periodic  orbits  around  collinear 
libration  points  (Fig.  1)  and  can  be  determined  using  a  dif¬ 
ferential  correction  algorithm  (Richardson,  1980).  The  ini¬ 
tial  conditions  of  these  halo  orbits  can  be  given  as 
x0  =  [x(0),y(0),z(0),i:(0),  y(0),z(0)]T.  Because  of  the 
orbit’s  symmetry  with  respect  to  the  (x  —  z)  plane,  we 
obtain  the  initial  conditions  y(0)  =  0,  x(0)  =  0  and  z(0) 
=  0.  A  differential  correction  algorithm  can  then  be  applied 
to  determine  the  other  initial  conditions.  The  detailed 
implementation  of  halo  orbits  can  be  found  in  Richardson 
(1980). 

Fig.  1(a)  and  (b)  illustrate  families  of  halo  orbits  around 
the  libration  point  L2  in  the  Sun-Earth  system.  The  ampli¬ 


1  , 

J=2jt  (ul  +  uy  +  uz)dt 


(12) 


where  ux,  uy  and  uz  are  the  control  variables  (acceleration) 
in  the  x,  y  and  z  directions,  respectively,  t0  is  the  initial 
flight  time  and  tf  is  the  duration  of  the  flight. 

By  introducing  the  new  variables  x\  =  x,  x2  —  =  z, 

x4  =  x,  x5  =  y  and  x6  =  z,  the  controlled  second-order  non¬ 
linear  dynamical  model  can  be  rewritten  in  the  state  space 
form  as  follows: 


Xi  =  x4 
X2  =  Xs 
X-i  =  x6 


(13) 

(14) 

(15) 


Fig.  1.  Families  of  halo  orbits  around  the  libration  point  L2  in  the  Sun-Earth  system:  (a)  northern  halo  orbits  and  (b)  southern  halo  orbits. 
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x5  =  -2x4  +  x2  -  (1  -  n)  -  ft  j2-  -  mM  t  ^  4 

rPS  rPE  rPM 


ax  ^3  -*3  -*3  . 

-/«)_ 3--^p - mM_^  +  uz 

rPS  rPE  rPM 


(16) 


(17) 

(18) 


where 

4  =  (*!  +fl)2  +  X12+X23  (19) 

IpE  =  (xi  —  1  +  jli)2  +  X^  +  Xj  (20) 

'pm  =  (*i  -  xm)2  +  (x2  -  y^)2  +  x^  (21 ) 

Meanwhile,  the  initial  and  final  conditions  of  the  opti¬ 
mal  transfer  trajectory  are  given  by 

x(to)  =  lXl(to),X2(to),X3(to),X4(to),Xs(to),X6(to)]  (22) 

x(tf)  =  [Mtf),X2(tf)>X3(tf),x4(tf),x5(tf),x6(tf)l  (23) 


The  initial  position  and  velocity  conditions  in  Eq.  (22) 
are  for  the  space  near  the  Earth.  The  final  position  and 
velocity  conditions  in  Eq.  (23)  are  for  the  desired  halo  orbit 
around  libration  point  L2  in  the  Sun-Earth  system. 

Thus,  the  optimal  control  problem  for  designing  the 
transfer  orbit  can  be  summarized  as  follows.  In  a  given 
fixed  time  interval  [<0,  f/],  a  spacecraft  with  the  dynamical 

model  given  by  Eqs.  (13)— (18)  is  controlled  from  the  initial 
conditions  Eq.  (22)  to  the  final  conditions  (23)  with  the 
minimum  energy  consumption  (12). 


3.2.  Generation  of  nominal  reference  trajectory  and  control 


With  the  nonlinear  optimal  control  model  established  in 
the  previous  section,  many  numerical  methods  (Betts,  1998; 
Conway,  2011;  Hull,  1997;  Biegler,  2010)  can  be  used  to 
solve  the  nonlinear  optimal  control  problem  and  generate 
the  reference  trajectories  and  control  histories.  These 
numerical  methods  are  divided  into  two  major  classes 
(Betts,  1998;  Conway,  2011):  indirect  methods  and  direct 


methods.  The  indirect  methods  transform  the  optimal  con¬ 
trol  problem  into  a  Hamiltonian  boundary  value  problem 
(HBVP)  by  calculus  of  variations  or  the  maximum  princi¬ 
ple;  the  HBYP  is  then  solved  using  numerical  methods, 
such  as  the  shooting  or  multiple-shooting  methods.  The 
primary  advantages  of  the  indirect  methods  are  the  high 
accuracy  of  the  solution  and  the  assurance  that  the  solution 
satisfies  the  first-order  optimality  conditions.  However,  the 
indirect  methods  have  several  disadvantages,  including 
small  radii  of  convergence,  the  need  to  analytically  derive 
the  HBVP  and  the  need  for  an  initial  guess  for  the  costate. 
In  the  direct  methods,  the  optimal  control  problem  is 
discretized  and  converted  to  large-scale  parameter 
optimization  problems  that  can  be  solved  using  nonlinear 
programming  methods  (Hull,  1997;  Biegler,  2010).  The 
direct  methods  are  advantageous  to  the  indirect  methods 
because  they  have  larger  radii  of  convergence  and  do  not 
require  an  initial  guess  for  the  costate. 

The  nonlinear  optimal  control  problem  has  the  essential 
characteristics  of  a  Hamiltonian  system.  The  fundamental 
property  of  Hamiltonian  systems  is  that  the  phase  flow  is 
a  symplectic  transformation.  Numerical  methods  that  pre¬ 
serve  the  symplectic  structure  are  more  effective  for  solving 
Hamiltonian  systems.  For  example,  the  symplectic  method 
exhibits  excellent  energy  behavior  and  accurately  reflects 
the  qualitative  behavior  of  the  solution.  Therefore,  we 
employ  a  recently  developed  symplectic  numerical  method 
(Peng  et  al.,  2011)  to  solve  the  nonlinear  optimal  control 
problem  in  this  paper.  The  nonlinear  optimal  control  prob¬ 
lem  is  transformed  to  nonlinear  algebraic  equations  by  the 
dual  variational  principle.  The  state,  control  and  costate 
variables  can  all  be  obtained  from  the  solutions  of  nonlin¬ 
ear  algebraic  equations.  The  detailed  realization  can  be 
found  in  Peng  et  al.  (2011). 

We  have  compared  several  strategies  for  the  design  of 
transfer  trajectories  to  test  the  effectiveness  of  the  symplec¬ 
tic  numerical  method  for  solving  the  optimal  control  prob¬ 
lem.  We  only  found  detailed  information  about  the  design 
of  transfer  trajectories  to  halo  orbits  around  the  Lx  point 
for  the  Sun-Earth  three-body  system,  so  we  compare  the 
results  of  our  method  with  the  dynamical  systems  method 
(Gomez  et  al.,  2005)  and  the  dynamical  system  method 


Table  1 

Parameters  for  generation  of  nominal  reference  trajectory  and  control. 


Parameter 

Value 

Units 

Mass  ratio  parameter  p 

3.034040E— 6 

Unitless 

Mass  ratio  parameter  mM 

3.694262E-8 

Unitless 

Orbital  parameter  RM 

3.942878E+5 

km 

Velocity  parameter  com 

1.547340E— 5 

rad/s 

Average  Sun-Earth  distance 

1.4959787E+8 

km 

Unit  time  length 

5.0225480E+6 

s 

Initial  position 

[1.4962200992,  0.0001404799,  0.0001228010]  x  108 

km 

Initial  velocity 

[4.4432029400,  2.3224844237,  -0.6322604756]  x  103 

m/s 

Final  position 

[1.5058570439,  0,  0.0056593377]  x  10s 

km 

Final  velocity 

[0,  4.588948873,  0]  x  102 

m/s 

Time  of  fight 

113-123 

day 
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Table  2 

Performances  of  the  transfer  trajectories  with  different  transfer  times. 


Transfer  time 
(days) 

:  Velocity  increment 

AFinitial  (m/s) 

Velocity  increment 
AFtransfer  (m/s) 

113 

14.5281 

133.4435 

115 

13.5297 

95.2373 

117 

12.9250 

71.6593 

119 

12.0538 

73.4220 

121 

11.5316 

112.3495 

123 

10.7767 

159.3128 

combined  with  the  optimal  control  method  (Serban  et  al., 
2002).  The  coordinates  of  the  departure  point  and  the  arri¬ 
val  point  are  taken  from  Tables  1  and  2  in  Gomez  et  al. 
(2005),  in  which  the  initial  condition  is  not  a  point  of  the 
stable  manifold  of  the  halo  orbit.  Using  these  parameters, 
zero  initial  velocity  error  and  transfer  time  173.25  days, 
the  simulation  was  performed  using  the  Serban  et  al. 
(2002),  Gomez  et  al.  (2005)  and  the  present  method.  The 
transfer  trajectory  and  nominal  halo  orbit  are  obtained 


using  the  present  method  and  are  plotted  in  Fig.  2.  The  fuel 
consumption  from  Gomez  et  al.  (2005)  is  14.1298  m/s,  the 
fuel  consumption  from  Serban  et  al.  (2002)  is  13.4831  m/s 
and  the  fuel  consumption  from  our  method  is  16.2790  m/ 
s.  The  fuel  consumption  from  our  method  is  slightly  higher 
than  in  the  other  two  methods,  which  may  be  a  result  of  the 
different  method  used  to  evaluate  fuel  consumption.  In 
contrast  to  the  instantaneous  velocity  changes  of  the 
dynamical  system  method,  the  control  inputs  of  our 
method  continuously  change.  Therefore,  the  evaluation  of 
fuel  consumption  in  our  method  may  have  some  differences 
from  the  other  methods.  Furthermore,  the  comparisons  of 
the  methods  demonstrate  the  effectiveness  of  the  symplectic 
numerical  method  for  solving  the  optimal  control  problem. 

The  nominal  reference  transfer  trajectory  with  the  initial 
conditions  given  in  present  Table  1  is  also  not  an  orbit  of 
the  stable  manifold  of  the  halo  orbit.  Fig.  3  shows  the 
transfer  trajectories  from  the  near  Earth  orbit  to  the  halo 
orbit  around  the  libration  point  L2  in  the  Sun-Earth 
system  obtained  using  the  parameters  in  Table  1  and  zero 


xio'3 


Fig.  2.  Transfer  trajectories  from  the  low  Earth  orbit  to  the  halo  orbit  around  the  libration  point  L\  in  the  Sun-Earth  system:  (a)  x-y  projection,  (b)  x-z 
projection,  (c)  y-z  projection  and  (d)  3D  space. 
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Fig.  5.  Hamiltonian  function  for  different  flight  ti 


initial  velocity  error.  The  different  solid  lines  plotted  in 
Fig.  3  show  the  different  nominal  reference  trajectories 
for  different  transfer  times,  i.e.  113, 115, ...,  123  days.  To 
fully  describe  the  performance  of  the  transfer  trajectory, 
the  initial  velocity  increment  A  Kinilia|  at  the  beginning  of 
the  trajectory  (near  the  Earth)  and  the  velocity  increment 
A  ^transfer  due  to  the  nominal  control  during  the  transfer 
trajectory  should  both  be  given.  However,  in  contrast  to 
the  instantaneous  velocity  changes  of  the  dynamical  system 
method,  the  control  inputs  of  our  method  change  continu¬ 
ously.  In  addition,  the  halo  orbit  transfer  missions  must 
generally  correct  the  launch  error  within  the  first  7  days 
after  launch  or  the  velocity  change  that  is  required  to  cor¬ 
rect  the  error  will  increase  beyond  the  spacecraft’s  capabil¬ 
ity  (Serban  et  al.,  2002).  Therefore,  we  evaluate  the  initial 
velocity  change  A  Fmitiai  within  3  days  after  launch  and 
the  velocity  change  AFtransfer  after  3  days  due  to  the  nom¬ 
inal  control.  The  detailed  performance  of  the  transfer  tra¬ 
jectories  with  different  transfer  times  are  given  in  Table  2. 

Furthermore,  the  feasibility  of  the  controlled  trajectory 
that  was  obtained  from  the  symplectic  numerical  method 
is  independently  validated  by  propagating  the  states  via  a 
separate  ordinary  differential  equation  solver  (the  ODE45 
function  in  Matlab).  The  relative  and  absolute  errors  are 
both  set  to  1  x  10-12  in  the  ODE45  function.  By  interpo¬ 
lating  the  values  of  the  control  variables  that  are  provided 
by  the  symplectic  numerical  method  and  then  integrating 
the  controlled  nonlinear  differential  Eqs.  (13)-(18),  the  rel¬ 
ative  deviations  of  the  controlled  trajectories  of  the  sym¬ 
plectic  numerical  method  are  compared  with  the  results 
of  the  numerical  integration  (Fig.  4).  The  results  show  that 
the  controlled  trajectories  obtained  using  the  symplectic 
numerical  method  are  similar  to  the  trajectories  from  the 
numerical  integration  method  and  demonstrate  that  the 
direct  results  differ  from  the  integration  results  by 


1  x  10-3  to  1  x  10-6.  To  confirm  the  validity  of  the  sym¬ 
plectic  numerical  method  for  solving  this  type  of  orbital 
transfer  problem,  different  Hamiltonian  functions  for  dif¬ 
ferent  times  of  flight  are  shown  in  Fig.  5. 

4.  Closed-loop  guidance  law  using  the  receding  horizon 
control  method 

The  purpose  of  designing  a  closed-loop  guidance  law  is 
to  correct  deviations  from  the  space  mission  trajectory  that 
are  produced  by  environmental  disturbances  and  other  per¬ 
turbations.  In  this  section,  the  closed-loop  guidance  law  is 
developed  based  on  the  receding  horizon  control  method 
and  the  known  reference  trajectory.  The  original  Riccati 
differential  equation  in  the  receding  horizon  control  is  then 
replaced  by  an  equivalent  convenient  form  of  the  Riccati 
differential  equation  that  is  based  on  the  generating  func¬ 
tion.  Finally,  a  high-efficiency  numerical  method  based 
on  the  generating  function  is  developed  to  solve  the  equiv¬ 
alent  Riccati  differential  equation  and  avoid  the  online 
integration  of  the  original  Riccati  differential  equation. 

4.1.  Closed-loop  guidance  law 

The  closed-loop  guidance  law  is  designed  based  on  the 
trajectory  error  with  respect  to  the  reference  trajectory 
and  the  receding  horizon  control  method.  Thus,  the  linear¬ 
ized  dynamical  model  given  by  Eqs.  ( 1 3)— ( 1 8)  with  respect 
to  the  reference  trajectory  is  derived  as 

<5x(?)  =  A(t)5x(t)  +  Bdu(t)  (24) 
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0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


(26) 


where  I3  and  03  are  the  three-dimensional  identity  matrix 
and  three-dimensional  zero  matrix,  respectively.  In 
addition, 
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Eq.  (24)  is  a  linear  time-varying  system  in  which  the  val¬ 
ues  of  the  coefficient  matrix  A  (7)  are  determined  by  the 
state  histories  of  the  reference  trajectory.  The  symbols  dx 
and  <5u  are  the  deviations  from  the  reference  trajectory 
and  control,  respectively.  Many  methods  can  be  used  for 
controlling  the  linear  time-varying  system  described  above; 
one  popular  approach  is  gain  scheduling.  In  this  method, 


controllers  are  designed  at  several  points  in  time.  The 
method  requires  the  repetitive  solving  of  algebraic  Riccati 
equations,  and  the  gains  are  then  interpolated  over  time. 
The  fundamental  issue  is  that  closed-loop  stability  cannot 
be  theoretically  guaranteed  by  such  a  gain-scheduled  con¬ 
troller  without  additional  conditions  (Lu,  2000).  Because 
the  linear  time-varying  system  is  continuously  regulated 
based  on  the  updated  information,  the  receding  horizon 
control  approach  is  more  robust.  Most  importantly,  the 
closed-loop  stability  has  been  proved  theoretically  (Lu, 
2000).  Therefore,  the  receding  horizon  control  approach 
is  used  to  control  the  linear  time-varying  system. 

The  receding  horizon  control  problem  at  any  fixed  time 
t  >  0  is  defined  to  be  an  optimal  control  problem  in  which 
the  performance  index 

/t-\-T 

[<5xt(t)Qc>x(t)  +  duT(T)Rdu(r)]dT  (27) 

is  minimized  for  a  chosen  time  interval  S  <  T  <  oo  (<5  is  a 
positive  constant)  subject  to  the  linear  time-varying 
dynamical  system  Eq.  (24)  and  the  constraint 
x(t+T)  =  0.  The  symbols  Q  e  IR',X"  and  R  €  IRmx'"are  the 
weighted  matrices.  The  weighted  matrix  Q  is  a  non-nega¬ 
tive  symmetric  matrix,  and  R  is  the  positive  definite  sym¬ 
metric  matrix. 

The  main  idea  of  receding  horizon  control  is  to  solve  the 
optimal  control  <5uoptm(f)  for  the  problem  described  above 
in  one  time  interval  [t,  t+T]  using  the  current  state  Sx(t) 
as  the  initial  condition.  Only  the  first  data  (5uoptm(f)  is 
employed  as  the  current  control  input  to  the  system;  the 
remainder  of  the  control  duoptm(f)  is  discarded.  For  the 
next  time  t,  the  solution  process  described  above  is 
repeated,  and  the  control  input  is  recomputed  (Lu,  2000). 
A  detailed  computation  of  a  control  input  obtained  using 
a  similar  receding  horizon  control  strategy  is  given  by  Lu 
(2000): 

<5u*(f)  =  — R~'BTP-1  (t,  t  +  T)Sx(t)  (28) 

where  P(f,  t  +  T)  satisfies  the  matrix  Riccati  differential  Eq. 
(29)  at  any  time  i  e  (t,  t  +  T)  =(t,  tT )  and 

P(t,  tr)  =  A(t)P(t,  tT)  +  P(t,  A)At(t) 

+  P(t,  fr)QP(T,  tT)  -  BR-‘Bt  (29) 

with  the  boundary  condition  P(fr,  tT)  =  P(f  +  T,  t  +  T) 
=  0. 

Assuming  that  the  matrix  pair  (A(f),  B(f)}  is  uniformly 
completely  controllable,  the  controlled  system  under  con¬ 
trol  law  (i.e.,  Eq.  (28))  is  uniformly  asymptotically  stable 
(Lu,  2000),  and  the  Riccati  differential  equation  has  a  sta¬ 
ble  solution.  From  a  practical  point  of  view,  the  online 
integration  of  the  Riccati  differential  equation  (i.e.,  Eq. 
(29))  for  every  time  t  7=  0  poses  a  significant  online  compu¬ 
tational  burden  if  the  control  strategy  Eq.  (28)  is  imple¬ 
mented.  The  key  problem  of  receding  horizon  control  is 
the  efficiency  of  real-time  computation.  Therefore,  in  the 
following  sections,  the  receding  horizon  control  is  devel- 
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oped  in  terms  of  the  generating  function  of  the  Hamilto¬ 
nian  system,  and  the  online  integration  of  the  Riccati  differ¬ 
ential  equation  is  avoided  to  increase  the  efficiency  in  real¬ 
time  computation. 

4.2.  Numerical  algorithm  for  solving  the  receding  horizon 
control  problem 

Park  et  al.  (2006)  stated  that  one  generating  function 
can  be  defined  for  the  linear  Hamiltonian  system  as 
follows: 


F„(f;t0) 

F/u(f;fo) 


Frffofo)  1  f<5xol 

-Fii(<;4>)J  l  8k  J 

(30) 


where  F xx(t;  t0),  F„(t;  t0)  =  F+r(f:  t0)  and  F u{t',  to)  are  coef¬ 
ficient  matrices  of  (F(x0,  A,  f;  t0),  8x0  is  the  initial  state  devi¬ 
ation  from  the  reference  trajectory  and  5k  is  the  costate 
variables  deviation.  The  generating  function 
(F(5x0,5k,t',to)  satisfies  the  following  Hamilton-Jacobi 
partial  differential  equation  (Park  et  al.,  2006): 


dF(5xQ,8k,P,t0) 

- of - +  H| 


«,(;»,)  A  =0 

d(5k)  J 


(31) 


and  the  initial  costate  5k&zamp-,0  and  state  variable  dxcan  be 
derived  from 

ax= -- Wp)  =  -F.((.(o)5xB+Fli((;4)i]  (32) 

Ao  jP(X°; \t;  t0)  F „(t;  t0)xoFx(t;  t0)5k  (33) 

0(xoJ 


where  H 

where  H  is  the  Hamiltonian  function  and  can  be  expressed 


Ij&wyT  q  atw  l/^wi 


H(5x,5k,t)=  =  <  /A 

V  ’  2  l<5A(f)J  A(t) 


I  m  J 


(34) 


According  to  the  Park  et  al.  (2006),  the  Hamiltonian 
function  (34)  can  be  formulated  with  the  generating  func¬ 
tion  using  Eqs.  (32)  and  (33)  as  follows: 

lTr— F^;i0)  01 

J  ij 

^rTia  «« 

By  substituting  the  generating  function  (30)  and  the 
Hamiltonian  function  (35)  into  the  Hamilton-Jacobi  par¬ 
tial  differential  Eq.  (31)  and  comparing  the  same  terms  of 
the  coefficient  matrices,  the  following  ordinal  differential 
equation  is  obtained: 

-  Fjtt(f;  to)  +  A(t)Fu(f,  to)  +  F u(f,  t0)AT(f) 

+  F^f;  f0)QF«(t;  t0)  -  BR  BT  =  0  (36) 


\  Q  AT(t)  1 
[a(0  -br  btJ 


Eqs.  (29)  and  (36)  have  the  same  form  and  are  both  Ricc¬ 
ati  differential  equations.  However,  Eq.  (36)  is  derived  dif¬ 
ferently  than  Eq.  (29)  and  is  expressed  with  the  generating 
function.  Therefore,  solving  Eq.  (36)  is  equivalent  to  solving 
Eq.  (29).  Once  the  solution  of  the  Riccati  differential  equa¬ 
tion  is  computed,  the  feedback  control  law  Eq.  (28)  of  the 
receding  horizon  control  problem  is  obtained.  The  solutions 
of  the  Riccati  differential  equation  (i.e.,  Eq.  (36))  are 
obtained  by  the  combination  and  separation  of  the  coeffi¬ 
cient  matrices  of  the  generating  function  {F(xo,  A,  t;  to)), 
which  replace  the  online  integration  process. 

In  view  of  Eqs.  (32)  and  (33),  the  relationship  of  the 
state  and  costate  can  be  expressed  in  terms  of  the  coefficient 
matrices  of  the  generating  function;  i.e.,  Eqs.  (32)  and  (33) 
have  an  equivalent  expression: 


m 


r  -F^fi;  to)  -  F u(t;  to)Fj  (t;  f0)F „(/;  t0) 

(<5xfi0n 

Wo)J 


-F u(f,to) 


(37) 


In  a  numerical  algorithm,  the  simulation  time  domain 
[to,  tf  \  is  divided  into  N  time  intervals  of  equal  time  step 
size  rj,  i.e.,  to  =  0,  t\  —  r\, . . .,  tk  =  kt], . . .,  t/=  Nt],  Hence, 
the  state  transfer  matrices  in  the  intervals  (tk-i,  tk),  {tk,  - 
tk+\)  and  {tk- 1;  tk+  i)  can  be  obtained  from  Eq.  (37)  as 
follows: 


r — (Fxc(4;  4— i|  +  Fa(4;  4-i)Fj/(4;  4-i)f«(4;  4-i ))  — Fju(4;  4-i)Fj  (4;  tk-i) 
[  -F-1  (4;  4-i)F«(4;  4-i)  -Fj  (4;  4-i) 


(38) 


[-(F,(4+i ;  4)  +  F;j.(4+i ;  tk)Fj  (4+1 ;  4)F»(4+i  ;  4))  — F;u(4+i  ;  4)Fr/  (4+1 ;  4) 
|_  — Fx/  (4+1 ;  4)F«(4+i  ;  4)  — F^1  (4+1 ;  4) 


fl>(4+i ;  4)  = 


(39) 
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x  [ — (F^(4+i  ;  4-i)  +  Fm(4+i  ;  4-iJF*1  (4+i ;  4-i)F„(4+i  ;  4-0) 


— F/u(4+i;  4-i  (4+i ;  4-i ) 


(40) 


According  to  the  properties  of  the  state  transfer 
matrices, 

®(4+i,  4-i)  =  0(4+1, 4)®(4, 4-i )  (41) 

Eqs.  (38)-(40)  are  substituted  into  Eq.  (41),  and  the 
combination  of  the  interval  coefficient  matrices  of  the  gen¬ 
erating  function  are  formulated  as 

Fxx(4+i;4-i)  =  Fxc(4;4_i) 

+  F^(4;  4-i)Wxi(4+i;  4)F^(4;  4-i)  (42) 

Fa(4+i;  4-i)  =  -Frf(4;  4-i)Wri(4+i;  4)  (43) 

Fju(4+i;4-i)  =Fii(4+i;4) 


every  time  step.  The  only  problem  for  the  solutions  of 
the  Riccati  differential  equation  is  the  unknown  coefficient 
matrices  of  the  generating  function.  Therefore,  the  rela¬ 
tionship  between  the  generating  function  and  the  state 
transition  matrix  is  derived,  and  the  key  point  of  solving 
the  generating  function  is  converted  to  the  solving  of  the 
state  transition  matrix. 

The  state  transfer  matrix  of  the  linear  Hamiltonian  sys¬ 
tem  (34)  can  be  given  by 


0(t;  t0)  = 


r®»(*;A>) 

L<M4*o) 


<M4*o)  1 

to)  \ 


(51) 


+  F^(4+i;  4)FU(4;  4-i)^F*t(4+i;  4) 

(44) 

Y  -  (I  +  F«(4+i ;  4)F«(4;  4-i))_1  (45) 

The  combination  formulas  (i.e.,  Eqs.  (42)-(44)  show 
that  if  the  coefficient  matrices  of  the  generating  function 
are  known  for  time  intervals  (4-4  4)  and  (4;  4+i),  then 
the  coefficient  matrices  of  the  generating  function  for  a  lar¬ 
ger  time  interval  (4_i;  4+1)  can  be  evaluated  by  the  combi¬ 
nation  formulas. 

Similarly,  according  to  the  properties  of  the  state  transi¬ 
tion  matrices, 

3>(4+i,4)  =  <&(4+i,4-i)«I>-1(4,4-i)  (46) 

Eqs.  (38)-(40)  are  substituted  into  Eq.  (46),  and  the  sep¬ 
aration  of  the  interval  coefficient  matrices  of  the  generating 
function  are  formulated  as 

F,jl(4+i;  4)  =  (^2Fm(4;  4_i)  —  Fx/l(4;  4-i))  ^^(4+1;  4-i) 

(47) 

Fxx(4+i;4)  =  —  Fxa(4+i;  4)Fx/(4+i;  4-i)^  (48) 

F^(4+i;4)  =  F/u(4+i;4-i) 

+ F^(4+i  ;  4-i  )F  (4 ;  4-i  )Fju(4;  4-i)F,*a  (4+1 ;  4) 

(49) 

n  =  (F«(4+i;4-i)  -  F„(4; 4-i))F^(4; 4-i)  (50) 

The  separation  formulas  (i.e.,  Eqs.  (47)-(49))  show  that 
if  the  coefficient  matrices  of  the  generating  function  are 
known  for  time  intervals  (4_i;  4)  and  (4-i;  4+i),  then 
the  coefficient  matrices  of  the  generating  function  for  a 
smaller  time  interval  (4;  4+1)  can  be  evaluated  by  the  sep¬ 
aration  formulas. 

Based  on  the  discussion  presented  above,  the  online  inte¬ 
gration  of  the  Riccati  differential  equation  is  replaced  by 
one  combination  formula  and  one  separation  formula  in 


which  relates  the  states  and  costates  at  time  t  with  the  ini¬ 
tial  states  and  costate  by 

r  <5x(f)  1  =  r  to)  <M4  to)  1  f  Sx(to)  ) 

l  5l(t]  J  L<M44)  «M4*o)J  1<M(4)  J 

The  state  transfer  matrix  of  the  linear  Hamiltonian  sys¬ 
tem  is  a  symplectic  matrix,  i.e. 

®T(t;4)J<D(t;4)  =  J  (53) 

where  J  is  a  unitary  symplectic  matrix. 

According  to  Eq.  (37),  the  state  transfer  matrix  can  be 
expressed  in  terms  of  the  generating  function.  Comparing 


Eq.  (37)  with  Eq.  (52)  gives 

Fxx(P,to)  =  <S>;{(f,to)<S>te(f,t0)  (54) 

Fx,(t;4)  -  -<&il(t:to)  (55) 

F«(i;  t0)  =  -®«(f;  t0)  +  ®**(4  fo)®u  (*;  to)®*  (#5  k)  (56) 
F u(t,  t0)  =  <M4  (t;  to)  (57) 


The  inverse  of  the  partitioned  matrix  (44)  exists  in 
Eqs.  (54)-(57).  The  length  of  the  two  adjacent  time  points 
(4  to)  cannot  be  made  arbitrarily  large  because  the  columns 
of  0;u(4  to)  become  more  linearly  dependent  as  the  interval 
(4  to)  increases  (Davison  and  Maki,  1973).  If  the  columns 
of  matrix  <D;j.(4  4)  are  linearly  dependent,  the  inverse  of 
matrix  <D^(4 10)  does  not  exist.  Consequently,  if  the  length 
of  interval  (4  4)  is  small,  the  inverse  of  matrix  <!>,+(  4 10) 
always  exists. 

Eqs.  (54)-(57)  relate  the  coefficient  matrices  of  the  gen¬ 
erating  functions  and  the  state  transfer  matrices  of  the  lin¬ 
ear  Hamiltonian  system.  Using  this  relationship  and 
based  on  the  symplecticity  of  the  state  transfer  matrix, 
we  employ  a  structure-preserving  numerical  algorithm 
(Hairer  et  al.,  2006)  to  solve  the  state  transfer  matrix 
and  then  to  solve  the  generating  function  accurately  and 
efficiently. 
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Fig.  6.  Simulation  process  of  a  spacecraft  with  a  guidance  law  transferring  to  a  halo  orbit. 


5.  Numerical  simulations 

In  this  section,  numerical  simulations  are  used  to  vali¬ 
date  the  effectiveness  of  the  guidance  law  that  is  based  on 
receding  horizon  control.  While  transferring  from  the 
near-Earth  orbit  to  the  halo  orbit,  the  spacecraft’s  devia¬ 
tions  are  corrected  by  this  guidance  law  under  the  assump¬ 
tion  of  a  bicircular  four-body  model.  The  detailed 
simulation  process  is  shown  in  Fig.  6. 

Fig.  6  shows  that  the  optimal  reference  trajectory  x*(t) 
and  the  reference  optimal  control  u*(t)  are  calculated  first. 
The  closed-loop  guidance  law  is  then  designed  based  on  the 
receding  horizon  control  and  the  linearized  dynamical 
model  from  Eqs.  ( 1 3)— ( 1 8)  with  respect  to  the  optimal 
reference  trajectory.  The  real  control  input  is  obtained  by 
combining  the  feedback  control  variable  c)u(t)  and  the 
reference  optimal  control  variable  u*(t);  i.e.,  u(t)  = 
u *(t)  +  Su(t).  The  real  state  x(f)  is  found  by  entering  the  real 
control  input  into  the  nonlinear  controlled  dynamical 
model  Eqs.  (13)— (18).  Finally,  the  deviation  Sx(t)  at  the 
next  time  point  is  generated  by  dx(t)  =  x(7)  -  x*(t).  This 
process  is  repeated  until  the  end  of  the  simulation. 

To  measure  the  fuel  expenditure  for  the  transfer  process, 
the  velocity  increment  AV  is  defined  as  (Kulkarni  et  al., 
2006) 

AVX  =  j‘f  \ux\dt,  AVy  =  J  '  \uy\At,  AVZ  =  f  '  \uz\dt  (58) 

AV=  AVX  +  AVy  +  AVZ  (59) 

where  the  fuel  expenditure  A  V  is  in  units  of  meters  per  sec¬ 
ond.  To  measure  the  final  relative  position  and  velocity  er¬ 
rors  with  the  closed-loop  guidance  law,  the  position  error 
Aep  and  the  velocity  error  Aev  are  defined  as 
Aep  =  max  (| \x(tf)  -  x*{tf)\,  \y{tf)  -  y\tf) |,  \z{tf)  -  z\tf)\) 

(60) 

Ae„  =  max(| x{tf)  -x*{tf) |,  \y{tf)  -|* (#/)[,  | z(tf)  -  z\tf)\) 

(61) 

where  x{tf),  y{tf)  and  z{tf)  represent  the  actual  final  posi¬ 
tion,  while  y*(tj )  and  z*(tj)  are  the  reference  final  po¬ 

sition  in  the  x,  y  and  z  directions.  The  final  velocity  error  is 
defined  similarly  in  Eq.  (61).  The  final  relative  position  er¬ 
ror  Aep  is  in  units  of  kilometers,  and  the  final  velocity  error 
is  in  units  of  meters  per  second. 


5.1.  Effect  of  initial  injection  error 

The  initial  injection  condition  is  important  for  the  preci¬ 
sion  and  quality  of  a  spacecraft’s  transfer  mission.  Two 
types  of  injection  error,  i.e.,  the  deterministic  injection 
error  and  the  stochastic  injection  error,  are  used  in  the 
numerical  simulations. 

For  the  deterministic  injection  error,  the  actual  depar¬ 
ture  velocity  from  the  orbit  near  the  Earth  has  been  mod¬ 
ified  using  the  following  formula  (Gomez  et  al.,  2005): 

,(0)=v;(1+ra)  (62) 

where  e  is  a  parameter  that  varies  between  0  m/s  and  10  m/ 
s,  v,)  is  the  reference  initial  departure  velocity  and  v(0)  is  the 
actual  departure  velocity. 

The  other  related  parameters  in  the  numerical  simula¬ 
tions  are  chosen  as  follows.  The  duration  of  the  flight  is 
120  days;  the  receding  horizon  control  parameter  is 
T=  12  days;  the  initial  position  is  given  in  Table  1.  The 
fuel  consumptions,  i.e.,  the  velocity  changes  obtained  by 
the  numerical  simulations,  are  given  in  Table  3.  AVX,  A  Vy, 
AVZ  and  A  V  are  velocity  changes  from  the  nominal  trajec¬ 
tory  due  to  the  nominal  control,  and  AVx,AVy,AVzandAV 
are  the  velocity  changes  due  to  additional  guidance  control 
for  the  initial  velocity  errors.  Table  3  shows  that  increases 
of  the  parameter  s,  i.e.,  the  errors  of  the  initial  velocity, 
increase  the  fuel  expenditures. 

The  transfer  trajectories  and  the  reference  optimal 
transfer  trajectories  obtained  by  increasing  the  initial 
velocity  error  to  100  m/s  and  increasing  the  receding  hori¬ 
zon  control  parameter  to  T  =  40  days  are  plotted  in  Fig.  7, 
and  the  actual  control  histories  and  the  reference  optimal 
control  histories  for  these  parameters  are  plotted  in 
Fig.  8,  in  which  the  dashed  lines  are  the  actual  transfer  tra¬ 
jectories  and  control  histories  with  initial  velocity  error, 
and  the  solid  lines  are  the  normal  reference  transfer  trajec¬ 
tories  and  control  histories  without  initial  velocity  error. 
The  actual  trajectories  diverged  from  the  reference 
trajectories  in  the  middle  portion  of  the  transfer  process. 
Therefore,  the  perturbation  of  the  initial  velocity  condition 
has  a  large  effect  on  the  transfer  process.  However,  the 
actual  trajectories  are  very  similar  to  the  reference  trajecto¬ 
ries  in  the  final  part  of  transfer,  which  indicates  that  a 
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Fig.  8.  Actual  control  histories  (dashed  line)  and  reference  optimal  control  histories  (solid  line). 
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Fig.  9.  Errors  from  the  optimal  reference  trajectory  based  on  stochastic  velocity  errors. 


tion  errors  between  —10  m/s  and  10  m/s.  One  hundred 
simulations  are  shown  in  Fig.  9.  All  of  the  deviations 
caused  by  the  stochastic  injection  errors  rapidly  converged 
to  zero.  Therefore,  the  guidance  law  proposed  in  this  paper 
based  on  the  receding  horizon  control  method  is  robust 
and  effective. 

5.2.  Effect  of  modeling  uncertainty 

In  the  previous  section,  the  present  guidance  law  was 
shown  to  be  robust  with  respect  to  initial  injection  errors. 
In  addition,  the  modeling  uncertainty,  which  includes  fac¬ 
tors  such  as  solar  pressure  and  other  ephemeris  dynamics, 
will  influence  the  orbit  transfer  mission.  Because  the  pres¬ 
ent  guidance  law  is  designed  based  on  a  linear  time-varying 


model  that  is  linearized  to  the  bicircular  four-body  model, 
the  complex  ephemeris  dynamics  model  will  generate  small 
errors  in  the  linear  time-varying  model.  To  validate  the 
robustness  of  the  present  guidance  law  with  such  uncer¬ 
tainty,  small  perturbations  from  the  reference  values  are 
added  to  the  time-varying  system  matrix  A (t)  in  Eq.  (24). 

The  flight  time  in  the  numerical  simulations  is  120  days, 
the  receding  horizon  control  parameter  is  T=  12  days  and 
the  initial  velocity  errors  are  all  10  m/s  (i.e.,  10  m/s  is  added 
to  the  initial  velocities  x(0),  j(0)  and  z(0)).  The  velocity 
changes  and  the  final  position  (velocity)  errors  with  the 
modeling  uncertainties  of  ±5%-±20%  obtained  from  the 
numerical  simulations  are  given  in  Table  4.  Figs.  10  and 
1 1  show  the  position  and  velocity  errors  from  the  optimal 
reference  trajectory  with  a  ±20%  model  uncertainty,  and 
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Table  4 

The  effect  of  model  uncertainty  for  velocity  increment  and  final  position/ velocity  errors. 


Model  uncertainty  AA 

Velocity  increment  AV  (m/s) 

Final  position  err 

■or  Aep  (km) 

Final  velocity  i 

0 

120.8811 

4.7228 

0.0025 

+5% 

116.8203 

3.5198 

0.0020 

-5% 

125.0188 

4.9349 

0.0026 

+10% 

112.8646 

5.8267 

0.0035 

-10% 

129.1957 

4.5905 

0.0025 

+15% 

109.1778 

3.1324 

0.0016 

-15% 

133.3755 

4.7528 

0.0025 

+20% 

110.1330 

5.0265 

0.0027 

-20% 

137.5282 

4.6498 

0.0025 

Fig.  10.  Position  errors  from  the  optimal  reference  trajectory  with  model  uncertainty. 


time  (days) 


Fig.  11.  Velocity  errors  from  the  optimal  reference  trajectory  with  model  uncertainty. 


2108 


H.  Peng  et  al.  I  Advances  in  Space  Research  51  (2013)  2093-2111 


Fig.  12.  Control  errors  from  the  optimal  reference  control  with  model  uncertainty. 


Table  5 

The  comparisons  between  the  proposed  method  and  the  backward  sweep  method. 


Parameter  T  (day) 

CPU  time  Tcm  (s) 

Velocity  increment  AV  (m/s) 

Final  position  error  Aep  (km) 

Final  velocity  ei 

rror  Aev  (m/s) 

GF 

BS 

GF 

BS 

GF 

BS 

GF 

BS 

1 

3.71 

41.47 

147.2161 

135.5798 

0.0328 

0.0291 

5.4475e— 5 

5.0922e-5 

2 

3.67 

41.65 

132.7274 

135.2843 

0.1203 

0.1385 

l.lOlOe— 4 

1.1673e-4 

3 

3.69 

41.85 

123.1537 

142.7726 

0.2127 

0.3038 

1.8595e— 4 

2.1881e-4 

4 

3.71 

41.65 

116.0734 

127.8883 

0.4590 

0.5649 

3.3765e— 4 

3.7627e-4 

5 

3.72 

45.31 

112.8173 

145.4822 

0.5847 

0.8615 

4.4841e— 4 

5.5053e-4 

10 

3.67 

46.20 

118.7869 

207.0493 

3.8201 

3.3551 

0.0021 

0.0019 

20 

3.78 

47.92 

126.8184 

162.7748 

12.3315 

11.1271 

0.0071 

0.0053 

30 

3.73 

53.73 

131.9519 

179.8914 

19.5670 

11.2277 

0.0175 

0.0212 

60 

3.74 

62.40 

140.6725 

197.0870 

298.1657 

713.5882 

0.1317 

0.4081 

120 

3.78 

99.98 

144.2283 

208.1440 

2018.40 

5705.01 

0.7529 

2.5981 

Fig.  12  shows  the  control  errors  from  the  reference  optimal 
control  with  a  ±20%  model  uncertainty. 

Table  4  shows  that  increasing  the  model  uncertainty, 
i.e.,  the  system  matrix  A(Y),  changes  the  velocities. 
However,  the  final  position  and  velocity  errors  are  all 
similar  to  the  normal  value.  Therefore,  the  present  guid¬ 
ance  law  can  overcome  the  impact  of  model  uncertainty 
and  satisfactorily  track  the  nominal  reference  trajectory 
and  control.  This  conclusion  can  also  be  obtained  from 
Figs.  ( 10)  ( 12);  with  a  ±20%  model  uncertainty  in  the 
system  matrix  A(t),  the  trajectory,  velocity  and  control 
input  errors  are  all  similar  to  the  nominal  reference  val¬ 
ues,  and  they  all  converge  to  the  nominal  reference 
values. 

5.3.  Effect  of  the  guidance  law  parameter 

In  the  two  previous  sections,  the  effectiveness  and 
robustness  of  the  proposed  guidance  law  were  validated 
using  numerical  simulations  that  included  initial  velocity 


errors  and  model  uncertainties.  The  efficiency  of  the  guid¬ 
ance  law  is  also  important  in  the  real  implementation  of 
transfer  missions.  Although  the  generating  function 
method  for  solving  the  receding  horizon  control  problem 
can  efficiently  solve  the  equivalent  form  of  the  Riccati  dif¬ 
ferential  equation  and  does  not  require  the  online  direct 
integration  of  the  original  Riccati  differential  equation, 
the  computational  efficiency  also  needs  to  be  validated. 
To  compare  the  ability  of  the  real-time  implementation 
of  the  receding  horizon  control  method,  a  Riccati  transfor¬ 
mation-based  backward  sweep  method  (Bryson  and  Ho, 
1975;  Ohtsuka  and  Fujii,  1997)  is  used  to  solve  the  same 
problem.  The  method  used  is  to  solve  the  Riccati  differen¬ 
tial  equation  at  every  time  step. 

The  parameters  in  the  numerical  simulations  are  chosen 
as  follows.  The  flight  time  is  120  days,  and  the  initial  veloc¬ 
ity  errors  are  all  10  m/s.  The  results  obtained  from  the  gen¬ 
erating  function  method  and  the  backward  sweep  method 
with  different  values  of  the  receding  horizon  control 
parameter  T  are  given  in  Table  5.  Table  5  also  shows  the 
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Fig.  13.  Errors  from  the  optimal  reference  trajectory  for  the  GF  (solid  line)  and  BS  (dashed  line)  methods. 
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Fig.  14.  Errors  from  the  optimal  reference  control  for  the  GF  (solid  line)  and  BS  (dashed  line)  methods. 


total  CPU  time  consumption  Tcpu,  the  velocity  change  A  V 
and  the  final  position  (velocity)  errors  Aep  (Ae„).  “GF”  and 
“BS”  are  abbreviations  for  “Generating  Function”  and 
“Backward  Sweep”,  respectively.  Many  conclusions  can 
be  obtained  from  Table  5.  First,  the  CPU  times  for  the 
GF  method  remain  nearly  constant  with  increases  of  T. 
However,  the  CPU  times  for  the  BS  method  increase  signif¬ 
icantly  with  increasing  values  of  T.  The  computational  effi¬ 
ciency  of  the  GF  method  is  at  least  one  order  of  magnitude 
better  than  the  BS  method  for  the  same  value  of  T.  Sec¬ 
ondly,  the  velocity  change  of  the  GF  method  decreased 
at  first  and  then  increased  with  increases  of  T.  The  velocity 
change  of  the  BS  method  has  no  obvious  statistical  charac¬ 


teristics  and  has  many  local  minima.  In  general,  the  GF 
method  gives  smaller  fuel  expenditures  than  the  BS  method 
for  the  same  value  of  T  (except  for  the  first  line).  The  GF 
and  BS  methods  give  final  position  (velocity)  errors  Aep 
(Aet,)  with  precisions  of  the  same  order  of  magnitude.  How¬ 
ever,  the  precision  of  the  GF  method  becomes  better  than 
the  BS  method  at  larger  values  of  T. 

Figs.  13  and  14  show  the  errors  from  the  reference  tra¬ 
jectory  and  control  for  T  =  10  for  the  GF  and  BS  methods, 
respectively.  The  solid  lines  denote  the  results  from  the  GF 
method,  and  the  dashed  lines  denote  the  results  from  the 
BS  method.  The  errors  from  the  optimal  reference  trajec¬ 
tory  and  control  for  the  GF  method  and  the  BS  method 
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converge  quickly.  However,  the  GF  method  has  smaller 
deviations  than  the  BS  method  in  the  convergence  process, 
especially  for  the  position  errors.  Finally,  the  maximum 
value  of  control  input  for  the  GF  method  is  smaller  than 
that  in  the  BS  method. 

6.  Conclusions 

This  study  examines  the  problem  of  transferring  a  space¬ 
craft  with  low  thrust  from  a  near-Earth  orbit  to  a  halo 
orbit  around  the  libration  point  L2  in  the  Sun-Earth  sys¬ 
tem.  A  guidance  law  based  on  the  receding  horizon  control 
method  is  proposed  to  study  the  guidance  problem  using 
the  nonlinear  bicircular  four-body  model.  Using  the  launch 
velocity  errors  and  the  uncertainties  of  the  nonlinear 
model,  the  guidance  law  is  implemented  to  compensate 
for  deviations  from  the  optimal  reference  trajectory.  The 
effect  of  the  initial  injection  condition  and  the  guidance 
law  parameters  are  studied  in  detail  using  numerical  simu¬ 
lations.  The  results  of  the  numerical  simulations  demon¬ 
strate  the  capabilities  of  the  proposed  guidance  law.  The 
generating  function  method  for  solving  the  receding  hori¬ 
zon  control  problem  has  an  advantage  over  the  backward 
sweep  method,  which  is  based  on  the  online  integration 
of  the  Riccati  differential  equation.  The  generating  func¬ 
tion  method  results  in  a  one  order  of  magnitude  improve¬ 
ment  in  computational  efficiency  over  the  backward 
sweep  method. 
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